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Abstract 

For a class of partial differential algebraic equations (PDAEs) of quasi-linear type 
which include nonlinear terms of convection type a possibility to determine a time 
and spatial index is considered. As a typical example we investigate an application 
from plasma physics. Especially we discuss the numerical solution of initial boundary 
value problems by means of a corresponding finite difference splitting procedure 
which is a modification of a well known fractional step method coupled with a 
matrix factorization. The convergence of the numerical solution towards the exact 
solution of the corresponding initial boundary value problem is investigated. Some 
results of a numerical solution of the plasma PDAE are given. 
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1 Introduction 



In this paper quasi-linear PDAEs for u = u(t, x), x e Q := (0, 1), of the form 

Au t + Bu xx + C[u]u x + Du =f(t, x), te (0, t e ), x G Q, (1) 

for some t e > are considered, u and / are mappings u, f : [0, t e ] x Cl — > 
M. n , n > 1, where / (supposed to be sufficiently smooth) is given. A, B, C[u] 
and D are real (n, n)— matrices where A, B and D are assumed to be constant. 
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All matrices may be singular, but A, B ^ 0. C[u] may depend on u. We 
suppose that it is linear in u. Typically, when C[u] is linear in u, vector C[u]u x 
describes physical convection. 

For system ([1]) we study classical initial boundary value problems (IBVPs). 
Initial values (IVs) may be decomposed 



u(0,x) = $ a (x) + $ c (x), x e fi (2) 
Uk{0,x) Uk(0,x) can be prescribed arbitrarily 
otherwise, 



Uk(0,x) Uk(0,x) cannot be prescribed arbitrarily 
otherwise, 



with (k = 1, n). The boundary values (BVs) are of similar form, 

u(t, x) =^f a (t,x) + #«.(*, x), te[0,t e ], xedtt = {0,l}. (3) 

This means the data which can be prescribed arbitrarily are in $ a , \l/ a . The 
consistent data (see e.g. [2] or [12]) are collected in <3> c , \l/ c . Furthermore, we 
assume that the compatibility relations 

$ a (x) + $ c (x) = * a (0, x) + ¥ c (0, x), x G 9n, (4) 

are satisfied. Because almost every PDAE has its own IV and BV distribution 
we avoid here to give a more detailed general description of these data. In 
section [5] we solve this problem for the plasma PDAE considered in the next 
section. 

This paper is organized as follows. In section [2] a nonlinear PDAE from physics 
is presented. In section [3] known general index concepts applicable also to 
nonlinear PDAEs are considered. The numerical solution of IBVPs by a finite 
difference splitting method is studied in section HI Especially, convergence 
results are given. A numerical example from plasma physics is presented in 
section [5j 



2 An application 



In this section we cite a mathematical model from plasma physics [3, Chap- 
ter 5] whose underlying system of equations is of type (jTJ). It describes the 
space-time-movement of a system consisting of ions (with mass density 
rii and positive electrical charge q) and electrons (with density n e and electrical 
charge —q) in the space domain Q. The charge q{n e — n») produces an elec- 
trical potential (f) such that the ions move under the corresponding electrical 
force —q4> x . The system of electrons is considered to be a gas with (constant) 
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temperature T e and pressure p = kBT e n e where ks is Boltzmann's constant. 
The relation between n e and <fi is given by the equilibrium of electrical and 
mechanical forces, 

qn e (j) x - k B T e n £)X = 0, 
and the equations of conservation of mass and momentum for the ions are 

( d d\ 

n ijt + (n#)i) x - Din iiXX = and m ; I — + Vi— I u< + q<p x = 0, 

respectively. Di is a (constant) diffusion coefficient, and Vi is the velocity of 
the ions. The equation for (p is Poisson's equation 

4>xx - 4:ir(n e - m) = 0. 

It is convenient to transform this system of equations by a linear transforma- 
tion into a new system of the form ([T]) with n = 4, / = and new dependent 
variables Ui and with matrices 

/ 1 \ / -6 \ / "2 ui \ /0000 

/i - oioo r — I oooo c\v] — I o ti 2 o di n_ oo o o 
A_ oooo j -°- o ooo ' b rJ ~ o o h ^ _ oooo 



-1 u 3 

oooo 



10-10 



(5) 



The new variables Ui are (up to a constant) given by U\ ~ rii, w 2 ~ Vi, u% ~ 
n e , U4 ~ 0. b > and o?i > are constants. 

Note that the matrices A, B are diagonal and singular, and C[u] is linear in 
the components of u (for short we say that C[u] is linear in u). 
Further examples of a system of type flTJ are a poroelastic model of a liv- 
ing bone [5H61TT0] and the incompressible Navier-Stokes equations if eq. ([1]) is 
generalized in obvious manner to two and three space dimensions. 



3 Indices 



First we introduce two definitions of classical linear spaces. Let I, m be non- 
negative integers, and let C^ m be the space of scalar real functions w(t,x), t G 
[0, £ e ], x G Q, whose derivatives up to the (I + m) th order (time derivatives 
up to the I th order and space derivatives up to the m th order) are continuous. 
With C l >™ the set 

C l 1 l rn :={w = (w 1 ,...,w ri ) T ; wteCj?, i = l,...,n] 

is defined. By C 1 ^ we denote a set of vector valued functions with vanishing 
BVs, 

Cjft := [w G ClT; w(t, 0) = w{t, 1) = o}. 



Published in Applied Numerical Mathematics 42 (2002) no. 1-3, pp. 297-314, 



doi: 10.1016/S0168-9274(01 )00157-X 



While the index of linear PDAEs has been considered by several authors, 
e.g. [2 3 4 8 l~2fT3] . the index for nonlinear systems is investigated little, see 
however [TT] where a discretization based index definition has been given, and 
|14j . We mention that we do not transform system (CO) to a first order system 
because we prefer in numerical calculations the original second order form. 
In this paper, a well known index concept for PDAEs is used to construct 
certain operators whose invertibility (if so) yields indices. The basic notion of 
a time index v t and a spatial index v x can be found, e.g., in [3|12fl4j . 
Throughout this paper the time index v t of ([T]) is of special interest. We 
assume that a solution u of the IB VP (JTJ) exists and that u is sufficiently 
differentiable. Furthermore, we suppose that u(t,x) = for x G dQ. When 
\l/ a and \l/ c (see eq. fl3])) are known, zero BVs can be obtained by a suitable 
transformation of u. 



3.1 Time index 

Definition 1 If the matrix A is regular, the time index v t of the PDAE (T7J) is 

defined to be zero. If A is singular, then v t is the smallest number of times the 
PDAE must be differentiated with respect to t in order to determine Ut as a 
continuous function oft, x, u and certain space derivatives of u— components. 

Since in most practical applications of PDAEs the time index is i/ t — 1 or 
v t = 2, we give here the formalism for these two indices for PDAEs of type 
([1]) only. First, an auxiliary result is stated. To simplify the notation we use 
in the following the summation convention (summation about twofold indices 
from 1 to n). 

Lemma 2 Let u be sufficiently smooth, and suppose that C[u] is linear in u. 
Then 

(I) d t C[u)u x = C^[u x ]u t + C[u]u tx where C\l\u x } := C ij)k u^ x , 
(II) d t CU[u x )u t = C^[u tx ]u t + CW[u x ]u tt . 



PROOF. Let Cij t k := -q^ 3 -, i,j,k G {1, Componentwise differentia- 

tion with respect to time yields 

8tCijUj tX = Cij t kUk,t u j,x + CijUjjx = (Cij t kUj,x) u k,t + CijUjjx- 

Statement (I) follows from the definition of . Using this definition again 
and the linearity of C in u we get 

u x] u k,t —Cij t k u j,tx u k,t + Cij^Uj,x u k,tt — C ik [Ut x ]Uk t t + C ik [UxlUhjt 
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which is (II) in component form. □ 



U2,x o o - 

u 2lX 

UA,x 

0- 



For example, C[u) given in fl5]) produces CWfuJ = 
With this Lemma and with definitions 

L :=Bd 2 x + D, M[u,u x ] := L + C (1) [u x ] + C[u]d. 



we find formally from eq. ([I]) under the assumptions of Lemma [2] after one 
time differentiation (this case is relevant for v t — 1) the derivative array 



(6) 



{ A °) 












v M[u,uJ A) 








Two time differentiations of 
/ 



JT]) produce the derivative array (for v t = 2) 



A 

M[u, u x ] 
2C^[u tx ] 





M[u, Ud 






.4 



V 



C[u]w, 







\ 



ft 

\ftt J 



(7) 



When A is singular, the coefficient matrices (denoted by A) of the systems 
©, dZD are singular in the sense that Av = where v ^ is an appropriate 
vector function (e.g. in the case of eq. ()6]), v = (uJ,uf t ) T , u G C% 2 )- However, 
the coefficient matrices might uniquely determine Ut- The analogous problem 
for linear time varying DAEs (differential algebraic equations) is discussed in 
[U p. 29], and for linear PDAEs of first order it is investigated in |14j . 
For the nonlinear PDAE ([T]), the derivative array or ([?]) is used to write 

Pu t =F (8) 

where P is according to eq. ([6]) or (J7J) an operator valued matrix, and the 
vector F does not depend on u t . First, we must find P and an appropriate 
domain of definition D(P). Second, the invertibility of P must be studied. 
We mention, that the matrix P in equation ([H]) is the analogue to the nonsin- 
gular diagonal matrix D(xj) defined in [T3], Definition 3.5. However, here we 
need not require a diagonal form of P. 



3.2 Spatial index 



If B is singular, we suppose that the PDAE is written in quasilinear form 
(with respect to second space derivatives), i.e. we transform B according to 
B = SqBSi 1 = ( f Q °) where So, Si are constant regular (n, n)— matrices. 
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Definition 3 The spatial index u x of a system with a regular matrix B is 
defined to be zero. When B is singular, v x is the smallest number of times the 
PDAE 

Au t + Bu xx + C[S7 l u]u x + Du = f(t, x) 

(u = S\u, A = SqASi 1 , and so on) must be differentiated with respect to x in 
order to obtain 

U ■ [U\ xx ^ U2 )XX i ••■) Um,xxi ^m+l,a;> ■••> ^"n,x) 

S . ' "> v ' 

m elements n—m elements 

as a continuous function of t, x, u, u t and ui jX , . . . ,u m;X . 

It is straightforward to derive for U a derivative array (by differentiations of 
the PDAE with respect to x) analogous to eq. © or (J7J). As a result, after a 
minimal number of x— differentiations one can find a representation of U of 
the form Q U = G where Q is an operator valued matrix. The vector G is 
independent of the components of U. This definition of the spatial index is 
according to the one given in [12]. It does not transform the PDAE (0Q) to a 
system of first order as is often done, e.g. in [14J. Such a transformation, in 
general, changes the indices. For example, one can show that the time-index 
here is if and only if the corresponding index in [H] is 1, a time- index 1 or 
2 here implies an index 2 there. 



3. 3 Time index of the plasma PDAE 



Here, we study the time index v t of the PDAE (TTJ under the assumption of zero 
BVs, u(t, 0) = u(t, 1) = 0, and / G Cf' 9 with suitable p, q. To be specific we 

1 2 

further assume that this IBVP has a solution u G C 4 ' which possibly exists 
only local in time. To determine the time index of system ([1]) we generate 
for ut - if possible - system (jSJ) with the condition that P defined on D(P) 
is invertible. In order to determine u t we differentiate the third and fourth 
equation of the PDAE with respect to t with the result 

-u 3 ,tx + U3,tU4,x + u 3 u 4 ,tx = and u 4 ,txx + u ljt - u S)t = 0. 

A new system can be formed by means of these two equations and the first 
two equations of the original system. This is a closed system of four equations 
for the components of u t in terms of u and derivatives which are not time 
derivatives of u— components. Obviously, the new system can be written 

-U2U2, x -d lU4iX ^ p 



1 


(1 


i) 


1 








[ 






1 


-1 


dl 



(9) 
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Now what is essential for the determination of the index of the PDAE is that 
u in P can be seen as a fixed element with the result that P : D(P) — > R(P), 
D{P) C CI 'o, is a linear operator (R(P) denotes the range of P). We can try 
to solve the first equation in (jUJ) for u t by standard linear theory. In particular, 
if P^ 1 exists we get v% as the least number of time differentiations which are 
necessary to get the equation for u% (in the example considered here, one 
differentiation with respect to t is needed). 

Remark 4 The right hand side of the first equation in (Tj|) does not play any 
role when v t is to be determined. This implies that v t is independent of f in 
equation (Qp. This is a reasonable result because the right hand side function 
should not influence the indices of a PDAE. 

We ask whether the linear operator P with D(P) C C 4 'o has for fixed u = 
u(t, x) an inverse. The answer comes from kernel N(P) whose elements z fulfill 

PZ — ( -Z3,x+U4,xZ3+U S Z 4tX I = 0. 



Z1-Z3+Z4 



xx 



This reduces to z\ = z 2 = and a linear homogeneous coupled system of two 
ordinary differential equations for Z3 and 2:4, —z^ x + u^ x z 3 + Usz^ x = and 
— z 3 + z^ xx = 0, with homogeneous BVs, -23(0) = and -2 4 (0) = z^l) = 0. 
z^ is also a solution of the equation — z^ xxx + u^ x z i>xx + ^32:4^ = which 
can be reduced to — y xx + u^ x y x + u 3 y = where y = y(t,x) := z± tX {t,x). 
By <fx and <£>2 we denote two fundamental solutions of the equation for y, i.e. 
y = K\ipi + K 2 ^2 is the general solution (Ki, % = 1,2, are independent of x). 
Then the following Lemma holds: 



Lemma 5 Let u be such that 



Then z 3 = z 4 = 0. 



V?m(0) ¥>2,x(0) 




^ 0. 



Since the proof is simple, it is omitted here (for details see [TO]). 

Under the assumptions of this Lemma we see that N(P) = {0}. Therefore, 

system can be solved for u t , and the time index is v t — 1. 

We mention that based on Definition [3] the spatial index of the plasma PDAE 

can be determined similarly. The result is v x = 0. 



4 Numerical solution by a finite difference method 



In this section we consider the numerical solution of IBVPs ([T]) - (j3J) by means 
of a fractional step difference method which is combined with a matrix fac- 
torization. The fractional step method and numerous variants of it are well 
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known, see e.g. [9]. By this method, the order of the system of equations which 
must be solved can be reduced considerably. The effort can be reduced further 
by a proper partition of the original system matrix (denoted by L below) into 
two splitting matrices (L = L\ + L 2 ). 

To describe the general procedure, we first rewrite the nonlinear part Cfti]^ 
of the PDAE (TjQ) as C[u, d x ]u which is more convenient sometimes. C[u, d x ] is 
an operator valued matrix. For example, with C[u] given in (jSJ) it is 



C[u]u x 



«2 Hi 

u 2 d x 
-1« 3 




U 2 ,x 
\U4,x ) 



d x u 2 
u 2 d x d\d x 
— d x u 3 d x 




\ ( Ul \ 

u 2 
u 3 

\U4j 



--: C[u,d x ]u. 



We suppose that A is singular and is given as A = l^ 1 ™ 1 °j where I ni is the 
unit matrix of order ri\ < n (ni > 1). Let n 2 —n — n\. Corresponding to this 
partition of A we introduce the notation u — u 2 ) T and M = (mII m]1 ) 
where the (n, n)— matrix M may be one of the matrices B, C, D. The compo- 
nent representation of Uk is u k = (w^, w^ nfc ) T , k = 1,2. According to this 
partition, PDAE (pQ) is written 









<B n 


B 12 \ 


/ 


:) 


u t + 




u x 




1° 






\B 2l 


B 22 J 


V 



C u [u,d x ] C 12 [u,d x ] 



u + 




Furthermore, let the operator valued matrices Lk[u], k = 1,2, be defined by 

(C kj = C k j[u,d x ]) 



( 





^2 i n i n \ ( u a2 



Ia[u] : 
L 2 [u\ : 

such that L[u] := Bd 2 x + C[u, d x ) + D = L^u] + L 2 [u 



\{B 21 d 2 £ + C 2l + D 2l ) {B 22 dl + C 22 + D 22 ) 

(B n d 2 x + C u + D u ) {B 2l 8l + C 21 + D 21 ) 




:i2i 



Remark 6 Sometimes, other operators L\, L 2 with L = Li + L 2 may be more 
convenient because the solution of the equations may be easier. In every case, 
the partition should be such that the identity AL 2 = L 2 does hold. This is 
needed in a factorization as explained below in eq. [To]) . 

In order to obtain an approximate numerical solution of IBVP (JTJ) - (j4j) we 
consider it with zero BVs by means of a difference method. The first time 
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derivative is approximated with an equidistant time step size r by 

u t (t m+1 ,x k ) « -(u m+1 (x k ) - u m (x k )), t m+1 = (m + l)r, m = 0, 1, ... 
r 

The corresponding IBVP is space discretized on an equidistant grid 

n h ■= { Xk = kh, k = 0, 1, M; h = 1/M, M > 1} 
by using difference formulas (k = 1, M — 1) 

5 2 1 

u xx (t,x k ) ~j^u k (t) := — (u k -i(t) -2u k (t) +u k +i(tyj, 

6 6 6 

u x (t,x k ) &—u k {t) or -j-u k (t) or —u k (t) (13) 

where 5q, 5 + and £_ are the usual central, forward and backward difference 
operators, respectively. Now suppose t G (0, t e ), t e > 0. We approximate eq. 
(JTJ) by the difference equation (m = 0, 1, k = 1, M — 1) 

m+l _ m 

A + L h [u k ]u k =f k . (14) 



Lhlu™} denotes a discretization of L[u(t m ,x k )\. We rewrite this equation as 

,m+l 

r 

and factorize approximately for r — > 



(A + tL.K]) fc - fc = - M«ZW (is) 



A + rL h [v%] « (A + rL ftl K])(J + rL M K]) (16) 
=4 + r(L w [<] + AL ft2 [<]) + 0(r 2 ) = A + rL fe [<] + 0(r 2 ), 



where L^i and L^ 2 are corresponding discretizations of (TTT]) and (TT2I) . respec- 
tively (we used AL h2 = L h2 and L h i + L h2 = L h ). u™ is considered to be an ap- 
proximation for u(t m , x k ). Therefore, we study for rn — 0,1, k — 1, M— 1 
fractional splitting 

(A + rL M K]) M r V2 =fr +1 -L h K]uT, (17) 

m+l _ m 

(I + rL h2 K]) Uk - % =nr +1/2 . (18) 
m° is (for every fc) given as IV. 

Remark 7 27m discussion shows that the fractional step method requires the 
solution of two linear systems of coupled equations, but each of them is, in 
numerous applications, of a considerably reduced order. 

The approximate factorization (TIB"!) implies the following Lemma. 
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Lemma 8 Suppose 



(1) v — v(t,x) G M n , t G (0,t e ), is for some t e > the exact solution 
(sufficiently smooth) of the IBVP §W\), 0-0, 

(2) the system of algebraic equations (T7\), / fi^j) has a unique solution. 



Then the method f77| ] ; U~R) and £/ie system l[T5\) are equivalent for r — > 0, i.e. 
eqs. [Tl\ ), {Tfy) approximate the original system to the same t— order as eq. 

m 



4-1 Convergence of the fractional step method 



We consider the convergence of the numerical solution calculated by the 
scheme (TH)) towards the exact solution v(t m ,Xk) of the corresponding IBVP 
for r — > and h — > under the condition (m+l)r = t (t fixed). The basic as- 
sumption is that C[v] = C^ + C 1 ^] is linear in v where the matrix C° = const. 
is chosen in such a manner that the matrices Go k , k = 1, ...,M — 1, defined 
below in fl39|) are regular. For example, a proper choice may be C° = C[v], 
C l [v] = C[v) — C[v] where vector v does not depend on t and x, e.g. v is a 
suitable mean value of v. 

First, we need the full truncation error defined for m — 0, 1, 

fc = l,...,M-lby 

m+1 _ m 

■m+1 /l fe As i r L.W|„.to+1 /m+1 /in"! 

a k :=A + L h [v k \v k - f k (19) 

where Lhlv™} := Bp + C[t>™]^ + D, and g = 1 (for one-sided differences) 
or q = 2 (for central differences). 5 stands for 5q or 5 + or <5_. Under the 
assumption that v is sufficiently smooth we Taylor expand v™, v™ ±1 in t m+ i 
and Xfc to get in lowest order with respect to r and h 

a™ +1 =A(vft 1 ~ \vT^r + 0(r 2 )) + B{v%£ + ^%£ji 2 + 0(/> 4 )) 

Since " wff-V] = C° + C\v™ +1 ] - rC 1 ^ 1 ] and 

Av^ 1 + Bv^ + (C° + C 1 ^ 1 ])^ 1 + = / fc m+ \ the truncation 

error can be written 

a™ +1 =AO rn+1 {T) + BO m+1 (h 2 ) + C[^ +1 ]O m+1 (/i g ) (20) 
+ C 1 [v k n + 1 ]O m+1 (rh q ) + C>™ +1 ]O m+1 (r)- 
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Here, e.g., O m+1 (r) = tw, w G M n , w independent of r. Second, scheme ([T 
is written by means of the Kronecker product in matrix form, 



i M -i®-+Qhiu m ])u 
Q h [u 



m+1 



A 



" ll : =—P®B 
at 



\p q ® C[U m ] + I M -i ® D 
qh 



m+1 



(21) 



where U r 



mT 

ur u 



mT \T 
1 j •••) "M-lJ > 



(/r r ,...,/s£ 1 ) T . and 



P := 



-2 1 
1 -2 1 



)(Af-l)x(M-l) 



1 -2 



(22) 



q and P g depend on the discretization of the first space derivative. For 
example, using central difference formula in ({TBI) , it is q = 2 and 



P :-- 



o 1 

-10 1 



>(M-l)x(M-l) 



-1 



(another differencing is chosen in fl35l) ). The (n(M — l),n(M — 1))— matrix 

P q <S> C[P m ] is a block matrix whose block at position (j,k) is P g jfc(7['uJ 1 ], 
u m £ R n Eqs> (pj) and pjj imply that ym+i satisfies equation 



.4 



.4 



Im-i ® ^+Qfc[^W m+1 = ® -JK m + P m+1 + >T +1 , 

^O^t) + E 2 O rn+1 (h 2 ) + E™ +l O m+1 {h q ) 
+ P 4 m+1 [O m+1 (r/i g ) + O m+1 (r)" 



(23) 
(24) 



where P x := I M ~i ® 4, P 2 := 7 M _! g> P, P^ +1 := 7 M _i g> C[V 3 ' +1 ], 

:= 7 M _! <g> and, e.g., O m+1 (r) means O m+1 (r) G R n ( M_1 ). 

The foregoing relations can be used to estimate a norm of the global error 
7] m := V m — U m . In the following we choose the discrete L 2 — norm defined for 

a vector v = (yi, f„(M-i)) T by := ZiEjfii *>k 
Subtracting eq. f l2"Tj) from eq. f l2"3"j) . we obtain 



4 

T 



Im-i ® - ) r/ m+1 + g ft [V m ]^ m+1 - Q h [tf m ]C/ m+1 = ( J M -i O 3 ) v m + A m+1 



7)- 



By the identity 
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the foregoing equation can be written 

G m V m+1 = (Im-i ® ^) V m ~ (Qh[V m ] - Q h [U m }) V m+1 + A rn+ \ (25) 

G m :=I M -i ®- + ^P®B + \p, ® C[U m ] + I M -i ® -D. (26) 
t hr qh 

Now we use the fact that 

Q h [V m ) - Q h [U m ] = \p q ® (C[V m ] - C[U m }) = \p q <g> C l [ V m ] 

qh qh 

(because the linearity of C[u] implies C[V m ] — C[£/ m ] = C 1 [r) m }). Furthermore, 
one can construct a (n(M — l),n(M — 1))— matrix C[V m+1 ] such that 

(P q g> C l \q m ]) V m+l = C[V m+1 }r] m . 

Therefore, eq. ([25]) takes the form {G~ m = (G m ) _1 ) 

v m+1 =G~ m \I M -i ® - - —C[V m+1 ]) 7] m + G- rn A m+1 
\ t qh ) 

or for short 

. = It follows for 771 = 0, 1, ... 

v m+i = #™#™-i...#<y + H m H m ~ l ...H 1 R 1 + ... + H m R m + R m+1 , 

\\v m+1 \\ <\\H m H m - 1 ...H°\\ \\r)°\\ + \\H m H rn - 1 ...H 1 \\ \\R l \\ (27) 
+ \\H m H m - l ...H 2 \\ \\R 2 \\ + ... + ||# m || \\R m \\ + \\R m+1 \\. 

Now we require for < mr = t G (0, i e ] that 



sup 

mgN 



{||/r n fr n - 1 ...ff J '||, j = l,...,m} < oo, (28) 



possibly under a restriction of one of the forms 



T T 
Ko < ~ < Kl, ^2 < T? < K 3- (29) 

ft AT 
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Assumption ( 128]) is discussed shortly in Remark [T2l With this condition, the 
right hand side in (|27|) can be estimated further to give 

m 

\\r] m+l \ \ < K Q \\rf\\ +K 1 V] \\R j+1 \ \ < K \\r}°\\ + {m + l)K x max \\R j+1 \\ 

j=0 je[o,m] 

< K \\r]°\\ +K 1 - max (30) 

r je[o,m] 

K , Kt, i?i are constants. Using \ \R j+1 \ \ = \\G- j A j+1 \\, relation (ED yields 

\\R j+1 \\ <\\G- 3 E^O^t) + \\G- ] E 2 \\0(h 2 ) + \\G- ] E{ +l \\0(h q ) (31) 
+ \\G-'Ei +1 \\[0(rh' 1 ) + 0(r)] 

or also 

||P i+1 || <||G?- J '||(||Si||0(r) + \\E 2 \\0(h 2 ) + ||^ +1 ||0(^) (32) 
+ \\Ei +1 \\[0(rh q ) + 0(r)]). 



Note that ineq. ( 13T|) is sharper than ineq. ( 132]) . Since the matrices E\ are block 
diagonal, their norms can be estimated easily. For example, E{ +1 is of the form 
Ei +1 = ^{CK +1 ],...,CK^ 1 ]}, and 

\ max {C T [v{ +1 ]C[v{ + %" 1 ' 2 

where X max (E) is the largest eigenvalue of E. Especially, if v is sufficiently 
smooth for some time t < t e , then all norms \\Ei\\ are bounded. Therefore, 
( 1321 can be written also 

\\R' +1 \ \ < K 3 \\G~ j \\ (0(t) + 0(h 2 ) + 0{h q ) + 0(rh q 



\\m +1 \\ = 


max 




k 



with some constant A3. 

An estimate of HC^ H can be based on the eigenvectors <pk of the symmetric 
matrix -j^P- By definition, -^Pcpk = ^k4>k, k — 1, M — 1. The eigenvalues 
are given by = — sin 2 (jjjj ■ They fulfill for h — > the asymptotic relation 
\ k = -(kir) 2 + 0(h 2 ) (because h = 1/M). The set {<p k} k = 1, M - 1} is 
assumed to be orthonormal. Let $ := yh{4>\, be the matrix of the 

eigenvectors and A := diag(Xi, Am-i)- Since $ _1 = $ T = $, it follows 
that (for short * n := $ <g) J n ) 

® B =tf n (A <g> fl)tt n , (33) 
/m-i ® + £>) =^ n (j M _! (8) OpA + D) ) * n . (34) 
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For simplicity, we assume that the first order derivative d x is discretized by a 
one sided difference approximation (q = 1), say 



h q ~h h 



-l 1 

-l i 



-l i 



= £ + H M - X ) . (35) 



Inserting this into eq. f l26|) and using eqs. (133]) and ( |34|) . G- 7 can be written 



G J =G + Gi, 



1 



/.w-i : ( ^ - ^G° + D ) A C 



Gi =- # M -i ® G° + P ® G 1 ^] 



(36) 
(37) 

(38) 



Eq. (|37|) implies that the (n(M — 1), n(M — 1))— matrix Go can be reduced to 
the set of (n, n)— matrices G fc, k — 1, M — 1, 



' Goi 



G = 



J*n, G oifc :=^-iG + P> + A fcJ B. (39) 



If the low order matrices Go/c, k — 1, M — 1, are regular (this is due to the 
choice of G°), then G is regular. Therefore, G J = G fin(M-i) + Gq X gQ , and 
it follows 



|G" J || < 



^n(M-l) + G GJ) || IIGq 1 ) 



(40) 



provided I n (M-i) + G 1 G{ is also invertible. The second factor on the right 
hand side can be written 



(Gq 1 !! = maxUGofc 1 !!,, 

k 



1 



max 



k [^min (Gg^Gofc)] ^ 



(41) 



where | \Q\ \ n is the spectral norm of a real (n, n)— matrix Q, and A m j n (GrQ k Gok 
is the smallest eigenvalue of the matrix G^Gofc. The first factor on the right 
of ineq. (HOI can be estimated by means of 



f-^n(M-l) + G 1 G{ 



< 



1 — || Gq G\ ' 
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provided HG^GiH < 1. It may be appropriate to estimate further 
HGo^iH < HGq 1 !! Il^ill. HGq 1 !! is given already in gU), and is 



\\G{\\ =-||-Jm_i ® C^W] + ® Cf^'] 



^^(II/m-i^c 1 ^': 



ft 
1 



Af-1 



max 




+ max 


C[u{] 




V fee[i,Af-i] 




n fce[l,M-2] 





(42) 



Here we used the identities 

((#m-i ® c[^']) T (# M -i ® c[w 



H M -i ® C[Lr- 



1/2 



A, 



V 



C M-2 C 'M-2 



1/2 



max 



1/2 



max ||C fe || n 



(for short = C Relations ( l4"Tj) and fT4"2"j) show that the norms 1 1 Cq -1 1 1 
and \ \G\\ \ can be estimated by norms of low order matrices (of order n where 
n is usually small, e.g. n = 4 for the plasma PDAE given in section [2]). 
We summarize the foregoing result in the following Lemma. 

Lemma 9 Suppose that for r G (0, r ] and h G (0, ho] (tq, h > 0), 
(m + l)r = t G (0,t e ], 

(1) the (ji,n) — matrices Gok, k = 1, M — 1, are regular, 

(2) \\u k \\ n < K , k — 1, M— 1, j = l,...,m + l, where K is some constant 
independent of r and h, 

(3) there is a positive constant 5q < 1 such that the condition 

WG^W \\G{\\ <( max \\G^\\ n ) \( max \\C l [ui}\\ n (43) 
" u " " " \fce[i,M-i] 11 Ufcl|n y ftVfce[i,M-i] v 7 



satisfied. 
Then\\G-i\\ < ^ 



+ max 

fce[i,M-2l 



G 1 



<5 



This condition plays a crucial role when proving the convergence of the dif- 
ference scheme considered here. We illustrate this by the following example. 



Example 10 Inequality is for the case n = 1, A = 1, B = —1, C = 
C° = const. < 0, D = and Dirichlet BVs of type of a CFL-condition known 
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from the discretization theory of hyperbolic differential equations. In this case 

T 

G = $ 



'M-l 



1 c° ~ 



1 C °\ r 
-I J M _! - A 



t h 

Therefore, Gofc = - — Hr — ^k which implies 



$, G x = —Hm-i. 
h 



mm 

k 



(l + |C°|^ + r|A fc | 



< T. 



Furthermore, H m _-J1m-\ ^ a diagonal matrix which has eigenvalues and 1, 
hence = 1, and 



\G 1 | I I \Gi\ \ < 



T- 



|C°| 



If this is required to be less than 1 (according to §4M), it resembles the well 
known CFL condition \C \r/h < 1. 

In order to prove convergence, the estimate ( }3Tj) should be applied to the error 
inequality ( 1301) . To estimate the norms | iG^'.E'zl |, / = 1, ...,4, we use again the 
decomposition (|36|) under the assumption that Gq 1 exists. Then the relation 

& = Go (/n(M-i) + G Q l G[) gives 



G-iE, 



< 



(jn(M-l) + G 1 G\^j G 1 E I 
1 



■n(M-l) 



Cj Lr 1 



G l E t 



Since E\ is (for all I) a block diagonal matrix, the representation (|37|) implies 
that also the norm ||G?~ 3 'i£||| can be estimated by the calculation of norms of 



(n, n)— matrices, 



\G- 3 Ei\\ < 



In(M-l) + G Q 1 G{ 



max HGofcll^max \\E,,y\\ n . 



Sometimes, this estimate may be useful. The result of the foregoing estimates 
is the following Lemma. 

Lemma 11 Suppose 

(1) C[u] = C° + C l [u] e W nxn is linear in u, 

(2) the exact solution v = v(t,x) of the IBVP is sufficiently smooth for 

\Wt\ \ < oo, 



t G (0,t e ], t e > 0, especially 

(3) (m+l)r = te (0, t e ], meN, t fixed, 

(4) \\U j \\ < oo, j = 1, ...m, 
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(5) 11^11 = 0, 

(6) su VmeN {\\H m H m ~ \..H% j = l,...,m} < oo, 

(7) inequality (f^ ) with 5 < 1 is valid. 

Then an upper bound of the global error (see ineq. (TffiTj) ) can be expressed in 
terms of norms containing Gq 1 only: 



jn+1 I 



K x t 
< — — < max 



1 — Sn r I fce[i,Af-l 



\G^A\\ n O(r) + \\G^B\\ n O(h 2 



j6[l,m+l] 



+ . nu« fl ( HGo^HOW + \\G^EI\\(0(t) + 0(r/i)) 



From this Lemma, one may get convergence results possibly under a restriction 

of the type given in ||G feA|| n , | (G^!!?! |„ and HGq 1 ^!!, i = 3,4, 

depend on r, ft and on the indices of the PDAE. For example, in the case 

/ o o o o \ 

of the plasma system [y t = 1), we have with the choice C° = o o -°i t) 1 

loo / 

||Gofe ly i||n = 0(r), \\GQ k B\\ n = 0(t), provided r and h 2 are related by 
r/h 2 < k where k is independent of r and h. By the way we note that in 
the plasma example ||Gofe IU is °f order 0{t 1 / 2 ) only. Therefore, one should 
not estimate ||Grj^4.|| n by the upper bound | \G$ k \ \ n \ \A\ \ n . The dependence 
of norms on the indices of the PDAE is discussed in [12]. 

Remark 12 



(1) The sixth requirement is a stability condition. Conditions of this type are 
well known in the theory of difference methods of time dependent partial 
differential equations (see also ffitf where the assumption is discussed for 
linear PDAEs without convection) . Unfortunately, a general easy method 
to verify the assumption when a convection term is present is not avai- 
lable. It should be studied separately for each problem of type (QP with 
given matrices A, B, C, D. 

(2) In most cases, the seventh assumption can be satisfied only when the step 
sizes r and h are restricted in a certain way and/or when the convection 
term is not dominant, see also example [7Z71 



5 Numerical example 



For an example we choose the plasma PDAE (with n — 4, d — 1) described 
in section [23 The matrices are defined in (jSj). Let the parameters in B and 
C[u] be defined by b = 0.02 and di = 1. First, according to (|2J) and (J3]) 
IVs and (Dirichlet) BVs must be specified. Since Ut = 1 (see section [373]), 
both terms $ a , $ c are different from zero. One can prescribe arbitrary IVs 
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for two components of u = u(t,x) only, say u 2 and M4 (these are the non 
zero components of $ ). One can show that we can assign especially the BV 
^3(0, 0) which is assumed to be non zero. Given this BV and 

g^{x) := u^(0,x) = K4 cos(27nr) 

(where the constant K4 is defined below), it follows from the third equation 
of the system (Tj[|), ([5]) that the consistent IV of u 3 is 

g 3 {x) :=u 3 (0,x) = u 3 (0,0)e 9 ^/e^ \ 

This can be inserted into the fourth equation of the system with the result 
that the consistent IV for u\ is 

9l (x) := Ul (0,x) = g 3 (x) - g A , xx {x) = u 3 (0, 0)e^ (x) /e^ 0) + An 2 g 4 (x). 

The IV for u 2 can be given arbitrarily, here we choose 

g 2 (x) := ^2(0, x) = K 2 x(x — 0.5), K 2 = const., 

and BVs are for t G It 

Ul (t, 0) = Ul (t, 1) = u 2 (t, 0) = 0, u 3 {t, 0) = u 3 (0, 0), u A (t, 0) = u 4 (t, 1) = K A . 

According to relation (TJ]), the BVs for u\ imply K 4 = — u 3 (0, 0)/(47r 2 ). With 
these relations the IBVP (JTJ) - (j3J) with matrices ([5]) is defined completely, 
and we may solve it by a finite difference method. The two parameters K 2 
and ^3(0,0) can be chosen arbitrarily, for an example let ^3(0,0) = 0.2 and 
K 2 = 0.4. Furthermore, let the two step sizes r and h be related by r = K h (or 
less restrictive by r < K h) where K is a positive constant (K = 0.5 below). 
This restriction has its origin in the second equation which is of hyperbolic 
type (inhomogeneous Burgers equation). In the example, the CFL condition 
for the second equation is satisfied for the data chosen because 

T 

-max{Kfc| mr= i} 

is much smaller than 1. Some values of this expression are given in Table 1 
in the line CFL 2 for different values of N = 1/h. Furthermore, in the table 

Table 1 



N 


20 


40 


80 


160 


320 


CFL 2 


0.0771 


0.0799 


0.0811 


0.0817 


0.0819 


ei 


0.0075 


0.0051 


0.0027 


0.0013 


0.0006 


e2 


0.0141 


0.0113 


0.0076 


0.0049 


0.0031 



ei is defined by 



U ith - U„ 



i,h/2\ 



1,2, where U^h and are the 
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numerical approximations of U{ at time t = 1 and at the space grid points for 
the two different step sizes h and h/2. Note that u\ and w 2 are the ion density 
and the ion velocity, respectively (see section |2j), and, therefore, u\ must be 
nonnegative. 

In Figure 5.1, a typical result of the finite difference solution of the IBVP with 
the data given at time t — 1 is shown. The components u\ and U2 which are 
the most interesting ones are presented. 

Fig. 1. 



6 Conclusion 



After an application of PDAEs with a nonlinear convection term we first consi- 
dered the determination of the time and spatial index of such systems. These 
were based on definitions given earlier in the literature, e.g. We reduced 
this problem to the question whether there is an inverse of some operator de- 
fined on a properly defined solution space. For the case of the plasma PDAE 
(a system with time index 1) this was shown in detail. In most practical ap- 
plications, PDAEs have time index less than 3. 

Then we studied the numerical solution of corresponding IBVPs by means 
of a finite difference splitting method familiar from the treatment of partial 
differential equations. Especially the convergence was considered for the case 
that both the time and space step sizes tend to zero. The main problem in 
these investigations is (in one space dimension) to handle the limit h — > in 
relation to the time step size r. For fixed h and r — > this problem is easy. 
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Some results concerning the problem when both h and r go to zero were given. 
Finally, some results of a numerical solution of an IBVP from plasma physics 
were presented. This example is interesting because the PDAE which is non- 
linear is a mixture of partial differential equations of parabolic, elliptic and 
hyperbolic type. 
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